Parametrization of coarse grained force fields for 
dynamic property of ethylene glycol oligomers/ water 

binary mixtures 

Tamio Yamazaki* 

Analysis technology center, Canon Inc. 3-30-20, Shimomaruko, Ota-ku, Tokyo 146-8501, 
&JQ[ Japan. 
^3 



^ • Abstract 

To evaluate shear viscosity of ethylene glycol oligomers (EGO) /water binary 
mixture by means of coarse-grained molecular dynamics (CG-MD) simu- 
lations, we proposed the self-diffusion-coefficient-based parameterization of 
non-bonded interactions among CG particles. Our parameterization proce- 
dure consists of three steps: 1) determination of bonded potentials, 2) scal- 
ing for time and solvent diffusivity , and 3) optimization of Lennard- Jones 
parameters to reproduce experimental self-diffusion coefficient and density 
data. With the determined parameters and the scaling relations, we evalu- 
ated shear viscosities of aqueous solutions of EGOs with various molecular 
weights and concentrations. Our simulation result are in close agreement 
with the experimental data. The largest simulation in this article corre- 
sponds to a 1.2 fis atomistic simulation for 100,000 atoms. Our CG model 
with the parameterization scheme for CG particles may be useful to study 
the dynamic properties of a liquid which contains relatively low molecular 
weight polymers or oligomers. 

Keywords: Coarse-grained molecular dynamics simulation, Ethylene glycol 
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1. Introduction 



Aqueous polymer solutions are widely used in industrial and household 
applications. For example, in ink-jet printing, polymer additives are used 
to control fixity and osmosis of the ink droplets to the paper. In order to 
estimate such properties, it is very important to understand the shear vis- 
cosity and the diffusivity of the polymer solutions. The molecular dynamics 
(MD) method at atomic level is the most common technique to estimate the 
shear viscosity of a solution or the self-diffusion coefficient of a component 
[H, [2I . However, even if a small oligomer of which the number of monomers 
is around 10, the longest relaxation time (for usual polymer solutions, it will 
be the rotational relaxation time of polymer chains), is several tens nanosec- 
onds, which is about 1000 times longer than the characteristic time of water 
molecules. The rotational relaxation time of the polymer chains increases in 
proportion to some power of its molecular weight. 

To obtain the ensemble averaged quantities, the sampling time (the time- 
average period) should be longer than the longest relaxation time of the 
system. The extreme long relaxation time complicates the estimations of the 
shear viscosity and the self-diffusion coefficient by using all-atom molecular 
dynamics (AA-MD) of which typical time-steps is femto-second order. 

In order to extend the accessible time scale, the coarse-grained (CG) tech- 
niques have been developed for several decades. The CG techniques, in 
which some atomic groups are represented by a single particle, are widely 
used for simulations of meso-scale phenomena in soft materials (lipids, sur- 
factants and polymers), for example, vesicle formation and fusion self- 
assembly (sl, H], micelle formation p, 0, M, 113] and pore formation in lipid 



bilayer |lll . Il2l . Il3l . Il4| . The coarse grained molecular dynamics (CG-MD) 
can achieve speed-up of up to several orders of magnitude faster than an 
AA-MD. Of course, the degree of speed-up depends on the details of the CG 
model or the system size. CG-MD method is very powerful tool to study 
of the static property, such as equilibrium structure of a large system, and 
it is also successful in investigating the fundamental mechanism of the long 
time-scale behaviors of soft materials. 

To apply this method to a quantitative estimation of the dynamics, it 
should be noted that the time in CG-MD trajectory is not equivalent to the 
time in the all-atom (AA) description, due to the lack of atomistic details in 
the CG model. An effective method to obtain the consistency of the molec- 
ular motions between AA-MD and CG-MD is to multiply the time scale of 



2 



the CG-MD by a constant factor. 

Recently, several authors reported the systematic studies of the dynamical 
and rheological properties of polymer systems polymer systems (polystyrene 
melt 15^3, polyamide-6,6 melt [16] and aqueous polyethylene glycol so- 
lutions 13 ) by the multiscale approach that combines AA-MD and CG- 
MD simulations through the scaling of time in the CG model. With the 
well-tuned potential energy function among the particles in the CG model, 
CG-MD gives us the reasonable possibility to investigate both statical and 
dynamical properties of large scale systems which include macromolecules. 
However, there are few studies about the parametrization for the non-bonded 
CG interactions, which can be adapted to the estimations of dynamic prop- 
erties of system, especially of the multi-component system. 

In this article, we will present the results of the CG-MD simulations, in- 
cluding the systematic determination of the non-bonded parameters for the 
ethylene glycol oligomer (EGO)/water binary mixtures based on the self- 
diffusion coefficient and density data. We will also present the results of 
the shear viscosity of the mixture estimated by means of the nonequilibrium 
CG-MD simulations with the newly determined force field parameters. 
The article is organized as follows. The experimental methods for measure- 
ment of the shear viscosity and self-diffusion coefficient are explained in the 
section two. The explanation of the CG model and the force-field parame- 
terization for EGO/water binary mixtures and the computational details are 
in the 3rd section. The results of the comparison between our simulations 
and the experimental measurements are shown in the Results and discussion 
section that follows. 



2. Experimental section 

2.1. Measurement of the shear viscosity and the self- diffusion coefficient 
2.1.1. Materials 

Reagent grade diethylene glycol (DEC), tetraethylene glycol (TEG) and 
PEG600 were purchased from Sigma Aldrich Chemical Company and used 
in this work without further purification. EGO/water binary mixtures were 
prepared gravimetrically using distilled water. D2O with a minimum deuter- 
ation degree of 99.95 % (Merck Co. & Inc., Darmstadt, Germany) was used 
for all experiments as the NMR lock solvent. 
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2.1.2. Experimental details 

The shear viscosities of the EGO aqueous solutions were measured at 293 
K by the RE80 viscometer manufactured by Toki Sangyo Co., Ltd. 
The self-diffusion coefficient measurements were performed at 293 K by pulse 
field gradient spin echo (PGSE)[l^ using standard ledbpgp2s sequence on 
the AvanceGOO NMR spectrometer (Bruker BioSpin GmbH, Rheinstetten, 
Germany). To avoid mixing H2O in the aqueous EGO solution and D2O 
as the NMR lock solvent, we used NMR tubes which have a double tube 
structure. In the NMR tube, the inner tube was filled with the NMR lock 
solvent D2O, and the outer tube was filled with aqueous EGO solution. 
In ^H-NMR spectra of aqueous EGO solutions, the peaks of around 5 ~ 3.6 
ppm are assigned to the protons of CH2CH2O groups of EGO, and the peaks 
of around 5 ~ 4.7 ppm are assigned to the protons of water molecule. The 
self-diffusion coefficients of assigned peaks are abbreviated in this article as 
Dego ~ 3.6 ppm) and Dqh ~ 4.7 ppm). Due to the proton exchange 
between hydroxyl groups of EGO end and water molecules, the diffusivity of 
EGO have a strong influence on the -Dqh- 

Under the assumption that the proton exchange between hydroxyl of EGO 
end and water molecules is very rapid, a self-diffusion coefficient of water 



(Dw) can be obtained by [19 



= -^Don - ^^Dego, (1) 

where Xi is the molar fraction of hydroxylic protons of EGO ends, and 
is the self-diffusion coefficient of water molecule. 



3. Theoretical section 

3.1. Coarse-graining of ethylene glycol oligomer and water 

In our coarse-grained model, EGO molecules and water molecules are 
represented by coarse-grained particles, as shown in Figure 1. 

EGO molecules are modeled by two types of particles ("PA") and ("PB"). 
The PA particle represents the oxyethylene unit of both ends of a ethylene 
glycol chain, and the PB particle represents the oxyethylene unit in a ethy- 
lene glycol chain. The mass of PA and PB particles are 53 amu and 44 amu 
respectively. A mass of an oxygen atom of ether is distributed to two ad- 
joined coarse grained particles at an equal ratio. 
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Figure 1: Coarse-graining scheme of a tetraethylene glycol and water. Water is modeled 
as a PW particle which corresponds to three water molecules. 

Water is modeled by a single particle (" P W" ) corresponding to three real wa- 
ter molecules. The mass of PW particle is 54 amu. DEG, TEG and PEG600 
molecules are modeled as PA-PA, PA-(PB)2-PA and PA-(PB)ii-PA, and we 
expressed these EGO species as EG02, EG04 and EG013 respectively. 

3.2. Coarse-grained pair potential 

We assume that the total potential energy for CG molecule is written as 

t/tot = + ^ang(eyfc) + ^non(i?.,), (2) 

where the terms, f/b, f^ang, f^non are the effective potential functions of bond 
length Ljj, bond angle Qijk-, and the distance between non-bonded CG parti- 
cles Rij respectively. Here, i, j and k are the indices of the CG particles. The 
non-bonded CG interactions Unon{Rij) are modeled using the 12-6 Lenard- 
Jones potential, 

where aij is the (finite) distance at which the inter-particle potential is zero, 
eij is the depth of the potential well. We use the Lorentz-Berthelot combi- 
nation rule for interaction between the different species, 

- r—. — - + (A') 

In the case the interaction between PB and PW particles, eij is calculated 
from the eq.(5) instead of eq.(4). 



e^j = 7 ■ V^i ■ ^3 



(5) 



3.3. Scaling relations 

As mentioned in Introduction, whenever the CG model is used, due to 
the lack of atomistic details, the speed of the time-evolution of the CG-MD 
trajectory is inconsistent with that of the AA-MD trajectory. To estimate 
of dynamical properties of the aqueous EGO solutions by means of CG-MD, 
we introduced the time mapping parameter S", which is defined by Kremer 



et.al. [15| as a ratio between the effective segment frictions in the AA model 
and in the CG model, 

^AA 

S = ^. (6) 

where C,^^ is an effective scalar friction coefficient of a segment in the AA 
model, and as C,^'^ is ones in the CG model. In this article, we have assumed 
that the ratio between the shear viscosities in the AA model rj^^ and in the 
CG model t]^'^ is equal to the parameter 5*, namely: 

^AA ^ ^ . ^CG^ (7) 

According to the Stokes-Einstein relation, the self-diffusion coefficient D is 
inversely proportional to r] and the hydrodynamic radius of a segment rh, 

D^{r^-r^)-\ (8) 

In our CG model, three water molecules are included in a single PW particle. 
The PW particle in the CG model has a larger hydrodynamic radius than 
the individual realistic water molecule, and thus the self-diffusion coefficient 
of water molecule in CG model is somewhat smaller than that in the AA 
model, even if we apply the time mapping parameter 5* properly. 
We define the hydrodynamic radius r^'^ of the molecule in CG model by 

r,c^ = ^ ■ r,^^, (9) 

where n is a number of atomistic molecules, which are included in one coarse- 
grained molecule, n = 1 (for EGO), n = 3 (for water) as shown in Figure 
1, and r^^ is the hydrodynamic radius of molecule in the AA model, the 
factor -^/n in eq.(9) comes from the assumption of the linear relationship 
between the cube of the hydrodynamic radius of a molecule and its volume. 
From eqs.(7) and (8), the scaling relation between self-diffusion coefficients 
of molecular systems, which are described by the AA model and by the CG 
model, is given by 
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where D^*^ is defined using the Einstein relation, 

D^^ = hm ^^{&) - rZliO)Y) , (11) 

where r^^{t) is coordinates vector of the center of mass of the molecule in 
CG model at time t, {■ ■ ■) denotes the ensemble average. In this article, S is 
dealt with as a constant value for simplification, though the large molecular 
weight dependence of S is argued in the CG-MD study of polystylene melt 
in Ref.15. 

In order to determine S, we consider the self-diffusion coefficient of pure 
water using CG-MD simulation -D^^er; ^^e experimental value -D^a^er 
observed by NMR measurement. From eq.(lO), S is given by 

S=^^-^. (12) 

water 

3.4- Parameterization of bonded potentials of EGO chain 

The Boltzmann inversion method is well-known as one of the techniques 



to evaluate effective mesoscale potential from atomistic simulation [20|. This 
technique enables us to determine the bonded intramolecular interaction 
especially. In this article, f/b(-^y) and f/ang(0ijfc) are determined by this 
method. In a thermal equilibrium system, we assume that an appearance 
probability P of a state vector of a system Q, obeys the Boltzmann distri- 
bution. 

PiQ) oc exp(-/3[/(Q)), (13) 

where U{Q) is a certain effective potential energy as a function of Q , and /3 
is l/ksT . 

Once P{Q) is obtained, a effective potential U{Q) can be straightforwardly 
determined from the inversion of eq.(13). In this article, we assumed that 
all of the bonded (stretching / bending) potentials in the EGO chain are 
approximately given by the same U, as we will show in eq.(15). 
To obtain P{Lij) and P{Qijk), the AA-MD simulation of a triethylene glycol 
dimethyl ether (TEGDE) in gas-phase is performed at 293 K. With our 
coarse-graining manner, a TEGDE molecule is modeled by three coarse- 
grained particles (PBp,p = 1, 2, 3) as shown in Figure 2. 

The center of PB particle is defined as the center of mass of the oxyethy- 
lene monomer unit, in which two oxygen atoms at both ends of the unit are 
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Figure 2: Coarse- graining scheme of a triethylene glycol dimethyl ether molecule. 



weighted by 0.5 and other atoms are weighted by 1.0. Histogram H{Lij) of 
bond length (bond PBi — PB2 and bond PB2 — PB3) and histogram H{Qijk) 
of bond angle(bond angle PBi — PB2 — PB3) are obtained from the 50-ns 
trajectory of the atomistic molecular dynamics simulation. Then, these his- 
tograms are renormalized by 

P(L.,) oc ^ , P(0.,.) cx (14) 

As a technical subject for the U{Q), because there are a lot of noises in 
the potential energy function obtained by the Boltzmann inversion scheme, 
it should be smooth by some appropriate functions. We assumed that the 
probability distribution function P(Q) can be expressed by the linear com- 
bination of Gauss functions Gi [2l|. Then the effective potential energy 
functions U{Lij) and U{Qijk) are obtained by 

m 

U{Q)/k^T = - ln(P(Q)) = - ln{J2 Gi) + const. 



1=1 

— m > ; exp 7^ + const. 



(15) 



where Ai, /i; and are total area, center position and width of the Gauss 
function Gi respectively, and m is the number of the Gauss functions for 
smoothing of P{Q). In this article, we decided m = 3 for the both potentials 
of bond length {Q = Lij) and bond angle {Q = Qijk)- The areas (Ai), center 
coordinates {fii) and widths {^i) are determined by curve fitting of the cal- 
culated bond length / bond angle distribution data. 
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3.5. Parameterization of Non-honded potentials 

In our CG model for the EGO / water binary system, there are three 
different particles (PA, PB, PW), and 7 parameters for the non-bonded in- 
teraction Unon{Rij), which are epw epA, ^pw, <^pw, o"pa, ctpb and 7. These 
non-bonded parameters for the Unon{Rij) are determined based on the ex- 
perimental data (the densities and the self-diffusion coefficients for several 
EGO / water binary systems). The parameterization procedure consists of 
the four steps listed below. 

STEP 1 : The parameter epw as a function of parameter apw is determined 
so as to reproduce the density of pure water at 293 K, which is 0.998 g/cm^ 



22| . In order to satisfy the experimental pure water density, epw is uniquely 



determined according to crpw The self-diffusion coefficient of pure water 



D^^^gj. is evaluated by CG-MD simulation with the fixed parameter values 
(epw and cxpw)- The time mapping parameter 5* is calculated from eq.(12). 
STEP 2 : The parameter epA as a function of parameter apA is determined 
so as to reproduce the density of pure diethylene glycol (EG02) liquid at 293 



K, which is 1.118 g/cm [23|. With the parameters for the PW determined 
in STEP 1, unique cxpA is obtained through minimizing of the error between 
the D'^^ calculated by CG-MD and the D^^ observed by NMR of the com- 
ponents in the EG02/water binary mixtures (EG02 weight fraction : 0.2, 
0.5 and 0.8). 

STEP 3 : For determination of the parameter epB, we used the density of 



pure tetraethylene glycol (EG04) liquid at 293 K, which is 1.125 g/cm^ 23 
and the self-diffusion coefficients of the components in the triethylene glycol 
(EG03)/water binary mixtures (EG03 weight fraction : 0.2, 0.5 and 0.8) 
observed by NMR. The EG04 and EG03 are modeled in this article as PA- 
PB-PB-PA and PA-PB-PA, respectively. With the parameters for the PW 
and the PA determined already in STEP 1 and STEP 2, the apB can be 
obtained by the same process in STEP 2. At the first execution of STEP 
3, 7 of eq.(5) is setted to 1. After performing STEP 4, 7 is revised to the 
optimized value. 

STEP 4: In the case of the calculating of the interaction between PB and 
PW particle, the parameter 7, which adjust the miscibility of the EGO chain 
in water, should be determined. The 7 is obtained through minimizing of 
the error between the D'^^ calculated by CG-MD and the D^^ observed by 
NMR of the components in the PEG600 (EG013)/water binary mixtures 
(EG013 weight fraction : 0.2). 
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STEP 3 and 4 are sequentially repeated until the proper ctpb and 7 are 
obtained. In STEP 2, 3 and 4, the root means of square errors (RMSE) are 
evaluated by 



RMSE 



1 ^ 

-J](logD,EXP_log^AA)- 



(16) 



where M is the number of data, is experimental self- diffusion coeffi- 

cients, D^^ is calculated self-diffusion coefficients from eq.(lO). 

3.6. Simulation details 

All simulations of this work are performed by GROMACS 4.0.5 ^]. In 
the atomistic molecular dynamics simulation of gas phase of TEGDE, the 



temperature is held at T = 293 K by Langevin thermostat j25| with cor- 
relation time r = 1.0 ps . The production time step for integration is dt 
= 1 fs. And cutoff radius for LJ and Coulomb potentials is 1.4 nm. The 
general Amber force field (GAFF) j26| is used as the atomistic force fields 
for TEGDE molecule. Atomic charges of a TEGDE molecule are assigned 



by AMl-BCC 271] method. These force fields and atomic charges are gen 



erated by antechamber 1.4 [28|, |26| and acpype 1.0 [29|. In the CG-MD 
simulations, to evaluate the density and the equihbrated liquid structure, 2 
ns MD simulations are performed in the NpT ensemble. Nose-Hoover ther- 



mostat [30|, |31|, l32] is used at 293 K to control temperature of the system. 
Par rinello- Rahman barostat 3^1 is used at 1 atm to control pressure of the 
system. In the 2 ns MD simulation, the instantaneous densities are calcu- 
lated from the last 1 ns trajectory every 1000 steps and then are averaged. 
Omitting the first 1 ns data of trajectory as relaxation time, the self-diffusion 
coefficient is calculated from 3 ns (for water, EG02/water, EG04/water bi- 
nary mixture) or 30 ns (for EG013/water binary mixture) MD simulation, 
which is performed at constant volume NVT ensemble, at the initial struc- 
ture is the last record of the former 2 ns MD simulation. 
Non-equilibrium molecular dynamics(NEMD) [i^] simulation is applied for 
the calculation of the shear viscosity of aqueous EGO solution. In order to 
calculate the shear viscosity, 20 ns (for water, EG02/water, EG04/water 
binary mixture) or 200 ns (for EG013/water binary mixture) MD simula- 
tion is performed. After dropped first 5 ns trajectory, the shear viscosity is 
calculated by analysis of NEMD trajectory. 
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There are 5 sets of MD/NEMD simulations with different initial structures 
and initial velocity profiles. The diffusion coefficients or the shear viscosi- 
ties are calculated from each MD/NEMD simulations, and then are aver- 
aged. The cutoff radius i?cut = 1-4 nm for non-bonded interaction among 
the coarse-grained particles and a production time step dt = 10 fs for inte- 
gration of Newton's equation are used as common conditions in all CG-MD 
simulations. 

4. Results and discussion 

4.1. Bonded interactions for PEG chain 

Parameters {Ai, fii and ^i) of eq.(15), which are determined from curve 
fitting of the data of U{Lij) and U{Qijk) are summarized in Table 1 and 
Table 2, respectively. 



Tabic 1: Parameters of bond length potential represented by eq.(15) 
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PA-PA,PA-PB and PB-PB 
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0.382 


0.023 


0.305 
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0.229 


0.020 


0.338 
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0.028 


0.018 


0.266 
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1)0U(1 au< 


!,!(' i)oU>utial re})reseiit(Hl In- ( 


-q.(f5) 


angle type 


1 


Ai 






PA-PB-PB and PB-PB-PB 


1 


238.840 


57.471 


190.567 




2 


45.375 


24.819 


123.986 




3 


31.826 


14.765 


101.560 



Figure 3 shows the comparison between the AA model and the CG model 
for the histograms of the length (a) and of the angle G^fe (b) of bonded 
PB particles of TEDME, as shown in Figure 2. 

4-2. Non-bonded interactions 

The Lennard- Jones parameters for PW, PA and PB particles, the pa- 
rameter 7 of eq. (5) and the time mapping parameter S are straightforwardly 
determined by using the procedures (STEP 1 to STEP 4) mentioned as sec- 
tion 3.5. In STEP 1, to satisfy the experimental liquid density of pure water 
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Figure 3: Probability distributions of bond length between bonded PB particles (a) and of 
bond angle among bonded PB1-PB2-PB3 particles (b) of TEDME. Open Symbols denote 
the distribution obtained by using AA-MD simulation. Broken line denotes the distribu- 
tion obtained by CG-MD simulation. 

(0.998 g/cm^) at 293 K and 1 atm pressure, epw is uniquely determined 
according to apw, which is hmited in range from 0.38 to 0.42 nm. If the 
parameter crpw is within this range, then the PW fluid has a stable liquid 
phase. Though the crpw can be selected arbitrary within this range, we se- 
lected (Tpw = 0.40 nm. The observed self-diffusion coefficient of pure water at 
273 K and 1 atm pressure is D^^^^ = 2.0 m^/s by using NMR. From eq.(12), 
5 = 6.19 is obtained from eq.(12). Through STEP 2 to STEP 4, the ctpa, 
(TpB and 7 parameters are one by one determined by minimizing the RMSE, 
which evaluated by eq.(16). These optimized non-bonded parameters are 
summarized in TABLE 3. 
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Table 3: Parameters of non-bonded potential represented by eqs.(3)-(5) 



apw 


O'PA 


crpB 


epw 




epB 


7 


0.40 


0.45 


0.46 


2.650 


4.356 


3.523 


1.13 



4-3. Diffusion coefficients and shear viscosities of EGO /Water binary sys- 
tems 

From eq.(ll), the self-diffusion coefficient is obtained by evaluating of 
the slope of the mean square displacement (MSD) of the center of mass of 
the EGO/water molecules. The MSDs were calculated for each EGO and 
water molecules in the simulation box using 3 ns (for water, aqueous EG02 
and EG04 solutions) or 30 ns (for aqueous EG013 solution) MD-simulation. 
Dropping the first 1 ns MSD data as relaxation time, the slope of MSD versus 
time was evaluated by linear regression. 

Figure 4a and 4b show the comparison between experimental and calculated 
self- diffusion coefficients. Figure 4a shows the self-diffusion coefficients of 
the EGO molecules plotted against the EGO weight fraction (W) of the 
EGO/water binary mixtures. Figure 4b shows the self-diffusion coefficients 
of the water molecules plotted in the same manner as Figure 4a. Three dif- 
ferent EGOs, which are EG02, EG04 and EG013, and three different EGO 
weight fractions (W= 0.2, 0.5 and 0.8) are presented, including the data of 
aqueous EG013 solutions (W = 0.5 and 0.8), which did not used in the 
parametrization of non-bonded potentials at Section 4.2. 
Experimental observations show that the self-diffusion coefficients (both of 
the EGO and the water) are linear on a log scale when plotted against W. 
The CG-MD simulation results are also linear as shown in Figure 4a and 
4b. The dependence of the D on W is increasing against the EGO molec- 
ular weight [Mw) as shown in Figure 4a, but the self-diffusion coefficients 
of the water molecule does not show such tendency (Figure 4b). The ex- 
perimental W and M^y-dependences of D are reproduced by the CG-MD 
simulations correctly. Additionally, we performed CG-MD simulation of 
EGO45(PEG2000)/water binary mixture {W = 0.2). We found that the 
calculated self-diffusion coefficients of EG045 and water are comparable to 
the experimental values based on NMR [l9j (data not shown), though these 
NMR measurements were observed at 298 K, differ from our CG-MD condi- 
tion T = 293 K. 

The shear viscosity was calculated by the nonequilibrium method described 
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Figure 4: Self-diffusion coefficients D of EGO (a) and water (b) molecules as a function 
of the weight fraction W of EGO in the EGO/water binary mixtures. The open circles, 
the open triangles and the open squares denote the D^-^^ values measured using NMR 
for EG02/water, EG04/water and EG013/water binary mixtures, respectively. The 
filled circles, the filled triangles and filled squares denote the D'~^'^ values calculated using 
CG-MD simulations for EG02/water, EG04/water and EG013/water binary mixtures, 
respectively. 



by B. Hess [3J], which estimates the shear viscosity of hquid from nonequi- 
hbrium simulation with an external cosine acceleration profile of the form 



k = — 



Acos {kz) 
27r 

I 



(17) 



where is the acceleration in the x direction, is the height of the simulation 
box, A is the magnitude of acceleration, z is the coordinate z-direction. The 
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shear viscosity of liquid r] can be estimated by 



2^ ^ ^ 
77 ^{t)=—^miVi,^{t)cos{kzi{t))/^mi, (18) 

P i=i i=i 

where p is the density of a system, rrii is the mass of the particle of index i, 
X is the velocity in the x direction of the particle of index i, N is the total 
number of particles in the simulation box, ri~^ is the reciprocal of viscosity 
of the system. The A parameter must be carefully selected, the shear rate 
should not be so high that the system is driven too far from equilibrium. 
The maximum share rate of the CG system sh^^ is 



sh^G -A-^— (19) 



For our simulations with: 1]'^'^ ^0.01 [kgm^^s^^], 1^ ^ 8 [nm], p ?«1000 
[kgm~^], and sh^^ is approximately 0.2 [ps nm~^] A. This shear rate should 
be smaller than one over the longest correlation time in the system. For usu- 
ally liquids, it will be the rotational correlation time of the largest molecule in 
the system. In the aqueous EG013 solution (W = 0.8), the rotational relax- 
ation time of end-to-end vector of the coarse-grained EGO 13 is approximately 
6000 ps. In this case, parameter A should be smaller than 0.001 [nm ps~^]. 
When the shear rate is too high, the observed shear viscosity will be too low. 
In this article, we used A — 0.0005 [nm ps~^] for all nonequilibrium CG-MD 
simulations. 1]^^ is estimated by eq.(9), with S — 6.19 and 77*-^*^ obtained 
from eq.(18). Figure 5 shows the comparison between experimental and cal- 
culated shear viscosities of the aqueous EGO solutions. For nine EGO/water 
binary mixtures evaluated in this article, the symbols are same as in Fig. 4, 
the calculated shear viscosities are agreed with the experiments. This means 
that the calculated values for two EGO-concentrations(iy = 0.2 and 0.5) 
are very close to the experimental data and overall data included W = 0.8 
shows the same tendency as the experimental results. In order to verify our 
CG model at larger molecular weight of EGO, we performed the additional 
calculations of the EGO22(PEG1000) /water and EGO45(PEG2000) /water 
binary mixtures. At the highest EGO-concentrations {W = 0.5 for EG022 
and 1^ = 0.4 for EG045) considered in this article, the rotational relaxation 
times of end-to-end vector of both the coarse-grained EG022 and EG045 
are lower than the ones of the coarse-grained EG013 in the aqueous EG013 
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solution {W = 0.8). Therefore, the parameter A in the eq.(18) and the sam- 
phng time length of ri~^ are same as the case of the EG013/water binary 
mixture {W = 0.8 ). Figure 6 shows the comparison between our calculated 
results and the experimental observations measured at 293 K [s^. We found 
that the calculated viscosities are comparable to the literature values for 
aqueous EG022 solutions. Although our calculated viscosities are slightly 
lower than the literature values for aqueous EG045 solutions, the tendency 
of the dependence of viscosity on the EG045 concentration seems to be in 
accordance with the experiments. 
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Figure 5: Shear viscosity 77 of the EGO and water binary mixture as a function of the 
EGO weight fraction W. The symbols are same as in Fig. 4. 
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Figure 6: Shear viscosity rj of the EGO and water binary mixture as a function of the EGO 
weight fraction W. The open circles and the open triangles denote the rj^^^ values in lit- 
erature [ail for EGO22(PEG1000)/water and EGO45(PEG2000)/water binary mixtures, 
respectively. The filled circles and filled triangles denote the ij'-^'^ values calculated using 
CG-MD simulations for EG022/watcr and EG045/water binary mixtures, respectively. 

5. Conclusions 

CG-MD simulations for the EGO/water binary mixtures were performed 
at 293 K and 1 atm pressure. The EGO chain was modeled by two types 
of particles, PA and PB, which represent the oxyethylene units of both ends 
and the middle of the chain, respectively. Also, three water molecules are 
included into a single PW particle. With our CG model, the number of par- 
ticles that should be considered in the simulation can be reduced ten-fold, 
and the time step for integration of Newton's equation increases ten times 
compared with atomistic simulations. 

The parameters for the CG model were determined by the systematic man- 
ner, as was shown in this article. The CG bonded potential parameters for 
the EGO chain were obtained by the Boltzmann inversion of the correspond- 
ing atomistic distribution functions. Due to the soft pair potential among 
the CG particles, the time-scale of CG-MD simulation is not equivalent to a 
realistic time scale. In order to estimate the proper dynamical properties by 
means of CG-MD simulations, the scaling relation for time in the CG model 
is introduced. Moreover, the hydrodynamic radius of PW particles in our 
CG model is larger than those of atomistic water molecules, due to the gain 
in volume of the PW particle as the result of the coarse graining of water. 
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as shown in Figure 1. Therefore, we also introduced the scahng relation for 
the water diffusivity based on the Stokes- Einstein law. 

With the bonded potential parameters and the scaling relations, the 12-6 
Lennard- Jones non-bonded potential parameters are straight-forwardly de- 
termined to reproduce the experimental observations (the density and the 
self-diffusion coefficient). We adopted the Lorentz-Berthelot mixing rule for 
the non-bonded interactions between the unlike particles. In this article, the 
Lorentz-Berthelot mixing rule is slightly modified, as shown in eq.(4) and 
eq.(5), to estimate the proper miscibility of the EG013 in water. By using 
the determined CG force-field parameters for the EGO/water binary mix- 
tures, our CG-MD simulation gives the estimations which agreed well with 
the experimental shear- viscosity data, including of the PEGlOOO/water and 
the PEG2000/water binary mixtures which were not used in our parame- 
terization procedure. The largest simulation in this article corresponds to 
a 1.2 /xs atomistic simulation for 100,000 atoms. Our CG model with the 
parameterization scheme for the CG particles may be useful to study of the 
dynamic properties of a liquid which contains relatively low molecular weight 
polymers or oligomers. 

In our future work, we plan to investigate the CG models for the sev- 
eral watersoluble polymers (e.g., polypropylene glycol, polyvinyl pyrrolidone, 
polyvinyl alcohol) , for the estimations of the shear- viscosity and the diffusiv- 
ity of these water mixtures through CG-MD simulations. 
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